
###
###
### Why are Subnational Governments Responsive?
###
###

library(rjags)

# a function for estimating 95% credible intervals for cutpoints
cutpoint <- function(x,y){
    mod <- glm(y~x,data=temp,family=binomial(link="logit"))
    inits <- data.frame(cp=coef(mod)[1]/coef(mod)[2],beta=coef(mod)[2])
    jdata <- na.omit(data.frame(x=x,y=y))
    jdata <- list(x=jdata$x,y=jdata$y,n=dim(jdata)[1])
    model <- jags.model(file="cutpoint.jags",inits=inits,
                       n.adapt=100,data=jdata,n.chains=1,quiet=T)
    output <- coda.samples(model, variable.names=c("cp"),
                     n.iter=1000,thin=1)[[1]]
    quantile(output[,1],prob=c(.025,.5,.975))
}

###
### Load the data
###

load(file="aux_data.RData")
load(file="indiv_data.RData")
st_est <- read.csv("states_estimates.csv")
st_est_pty <- read.csv("states_parties_estimates.csv")

# drop DC and missing states
data <- data[! is.na(data$state) &  data$state != "District of Columbia",]

###
### Estimate Indifference Points by Year
###

states <- unique(data$state)
years <- unique(data$year)
st_data_year <- as.data.frame(matrix(NA,length(years)*length(states),24))
names(st_data_year) <- c("state","year","SEN","SEN025","SEN975","GOV","GOV025","GOV975","AG","AG025","AG975","SEC","SEC025","SEC975","SA","SA025","SA975","SS","SS025","SS975","DEM","DEM025","DEM975","N")

minN <- 20

for (i in 1:length(states)){
    print(states[i])
    for(j in 1:length(years)){
        print(years[j])
        ind <- (i-1)*length(years) + j
        st_data_year$state[ind] <- states[i]
        st_data_year$year[ind] <- years[j]
        temp <- data[data$state==states[i] &
                     data$year==years[j],]
        if(length(which(! is.na(temp$SEN))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SEN)
        st_data_year$SEN025[ind] <- cps[1]
        st_data_year$SEN[ind] <- cps[2]
        st_data_year$SEN975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$GOV))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$GOV)
        st_data_year$GOV025[ind] <- cps[1]
        st_data_year$GOV[ind] <- cps[2]
        st_data_year$GOV975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$AG))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$AG)
        st_data_year$AG025[ind] <- cps[1]
        st_data_year$AG[ind] <- cps[2]
        st_data_year$AG975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$SEC))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SEC)
        st_data_year$SEC025[ind] <- cps[1]
        st_data_year$SEC[ind] <- cps[2]
        st_data_year$SEC975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$SEC))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SA)
        st_data_year$SA025[ind] <- cps[1]
        st_data_year$SA[ind] <- cps[2]
        st_data_year$SA975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$SEC))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SS)
        st_data_year$SS025[ind] <- cps[1]
        st_data_year$SS[ind] <- cps[2]
        st_data_year$SS975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$DEM))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$DEM)
        st_data_year$DEM025[ind] <- cps[1]
        st_data_year$DEM[ind] <- cps[2]
        st_data_year$DEM975[ind] <- cps[3]
        }
        st_data_year$N[ind] <- dim(temp)[1]
    }
}

# zero out the few cutlines where 
# precision is very low
st_data_year$GOV[st_data_year$GOV975-st_data_year$GOV025 >1] <- NA
st_data_year$SEN[st_data_year$SEN975-st_data_year$SEN025 >1] <- NA
st_data_year$AG[st_data_year$AG975-st_data_year$AG025 >1] <- NA
st_data_year$SEC[st_data_year$SEC975-st_data_year$SEC025 >1] <- NA
st_data_year$SS[st_data_year$SS975-st_data_year$SS025 >1] <- NA
st_data_year$SA[st_data_year$SA975-st_data_year$SA025 >1] <- NA
st_data_year$DEM[st_data_year$DEM975-st_data_year$DEM025 >1] <- NA

# merge estimates to estimates of state-level opinion
st_data_year <- cbind(st_data_year,st_est[match(st_data_year$state,st_est$name),])
st_data_year <- cbind(st_data_year,st_est_pty[match(st_data_year$state,st_est_pty$name),])

# save the results just in case
save("st_data_year",file="state_data.RData")
# load(file="state_data.RData")

##
## Estimate Pooled Indifference Points
##

st_data <- as.data.frame(matrix(NA,length(states),24))
names(st_data) <- c("state","year","SEN","SEN025","SEN975","GOV","GOV025","GOV975","AG","AG025","AG975","SEC","SEC025","SEC975","SA","SA025","SA975","SS","SS025","SS975","DEM","DEM025","DEM975","N")

minN <- 20

for (i in 1:length(states)){
    print(states[i])
        ind <- i
        st_data$state[ind] <- states[i]
        temp <- data[data$state==states[i],]
        if(length(which(! is.na(temp$SEN))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SEN)
        st_data$SEN025[ind] <- cps[1]
        st_data$SEN[ind] <- cps[2]
        st_data$SEN975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$GOV))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$GOV)
        st_data$GOV025[ind] <- cps[1]
        st_data$GOV[ind] <- cps[2]
        st_data$GOV975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$AG))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$AG)
        st_data$AG025[ind] <- cps[1]
        st_data$AG[ind] <- cps[2]
        st_data$AG975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$SEC))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SEC)
        st_data$SEC025[ind] <- cps[1]
        st_data$SEC[ind] <- cps[2]
        st_data$SEC975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$SEC))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SA)
        st_data$SA025[ind] <- cps[1]
        st_data$SA[ind] <- cps[2]
        st_data$SA975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$SEC))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$SS)
        st_data$SS025[ind] <- cps[1]
        st_data$SS[ind] <- cps[2]
        st_data$SS975[ind] <- cps[3]
        }
        if(length(which(! is.na(temp$DEM))) > minN){
        cps <- cutpoint(x=temp$ideal_point,y=temp$DEM)
        st_data$DEM025[ind] <- cps[1]
        st_data$DEM[ind] <- cps[2]
        st_data$DEM975[ind] <- cps[3]
        }
        st_data$N[ind] <- dim(temp)[1]
}

# zero out the few cutlines where 
# precision is very low
st_data$GOV[st_data$GOV975-st_data$GOV025 >1] <- NA
st_data$SEN[st_data$SEN975-st_data$SEN025 >1] <- NA
st_data$AG[st_data$AG975-st_data$AG025 >1] <- NA
st_data$SEC[st_data$SEC975-st_data$SEC025 >1] <- NA
st_data$SS[st_data$SS975-st_data$SS025 >1] <- NA
st_data$SA[st_data$SA975-st_data$SA025 >1] <- NA
st_data$DEM[st_data$DEM975-st_data$DEM025 >1] <- NA

# merge estimates to estimates of state-level opinion
st_data <- cbind(st_data,st_est[match(st_data$state,st_est$name),])
st_data <- cbind(st_data,st_est_pty[match(st_data$state,st_est_pty$name),])

# merge Caughey Warshaw and Shor McCarty data
st_data <- st_data[! is.na(st_data$state) &  st_data$state != "District of Columbia",]
st_data <- cbind(st_data,aux_data)

# save the results just in case
save("st_data",file="state_data_pooled.RData")
# load(file="state_data_pooled.RData")

######
###### Analysis
######

###
### Figure 1: state policy v state opinion
###

st_data$policyeconlib_est <- st_data$policyeconlib_est*sign(cor(st_data$policyeconlib_est,
                                                                st_data$raw_mean))
pdf("basic_representation.pdf")
plot(policyeconlib_est~raw_mean,data=st_data,axes=F,
     xlab="State Ideal Point",ylab="State Economic Policy")
axis(1)
axis(2)
abline(lm(policyeconlib_est~raw_mean,data=st_data))
dev.off()

###
### Figure 2: State Partisan Indifference Points (pooled over years)
### 

st_data <- st_data[! is.na(st_data$state) & st_data$state != "District of Columbia",]

pdf("state_parties_all.pdf",width=7,height=10)
par(mar=c(5.1,10,2.1,2.1))
plot(NA,ylim=c(0,50),xlim=c(-.6,.6),axes=F,
     ylab="",xlab="party indifference point")
abline(h=c(1:50),col="grey",lwd=.5)
axis(1)
states <- unique(st_data$state[order(st_data$raw_mean)])
axis(2,labels=states,at=1:length(states),las=1)
for(i in 1:length(states)){
    state <- states[i]
    print(state)
    ind <- which(st_data$state==state)
    xs <- st_data$DEM[ind]
    ys <- rep(i,length(xs))
    segments(st_data$DEM025[ind],ys,st_data$DEM975[ind],ys,lwd=1)
    segments(xs,ys-.4,xs,ys+.4,lwd=2)
    points(st_data$raw_mean[ind][1],i,col="black")
}
abline(v=median(st_data$DEM,na.rm=T),lty=2)
dev.off()


###
### Figure 3: the relationship between indifference points for various offices
### and state idealpoints
###
par(mfrow=c(1,1))
pdf("offices_idealpoint.pdf",width=7,height=7)
plot(st_data_year$GOV~st_data_year$raw_mean,xlim=c(-.8,.8),ylim=c(-.8,.8),
     axes=F,xlab="state idealpoint",ylab="governor indifference point")
segments(st_data_year$raw_mean,st_data_year$GOV975,st_data_year$raw_mean,st_data_year$GOV025,lwd=.5)
abline(lm(st_data_year$GOV~st_data_year$raw_mean))
axis(1,at=seq(-.8,.8,.2));axis(2,at=seq(-.8,.8,.2))
plot(st_data_year$SEN~st_data_year$raw_mean,xlim=c(-.8,.8),ylim=c(-.8,.8),
     axes=F,xlab="state idealpoint",ylab="senator indifference point")
segments(st_data_year$raw_mean,st_data_year$SEN975,st_data_year$raw_mean,st_data_year$SEN025,lwd=.5)
abline(lm(st_data_year$SEN~st_data_year$raw_mean))
axis(1,at=seq(-.8,.8,.2));axis(2,at=seq(-.8,.8,.2))
plot(st_data_year$SA~st_data_year$raw_mean,xlim=c(-.8,.8),ylim=c(-.8,.8),
     axes=F,xlab="state idealpoint",ylab="state assembly indifference point")
segments(st_data_year$raw_mean,st_data_year$SA975,st_data_year$raw_mean,st_data_year$SA025,lwd=.5)
abline(lm(st_data_year$SA~st_data_year$raw_mean))
axis(1,at=seq(-.8,.8,.2));axis(2,at=seq(-.8,.8,.2))
plot(st_data_year$SS~st_data_year$raw_mean,xlim=c(-.8,.8),ylim=c(-.8,.8),
     axes=F,xlab="state idealpoint",ylab="state senate indifference point")
segments(st_data_year$raw_mean,st_data_year$SS975,st_data_year$raw_mean,st_data_year$SS025,lwd=.5)
abline(lm(st_data_year$SS~st_data_year$raw_mean))
axis(1,at=seq(-.8,.8,.2));axis(2,at=seq(-.8,.8,.2))
dev.off()

###
### Figure 4: party midpoint v. indifference points
###

st_data$hou_cl <- (st_data$hou_dem + st_data$hou_rep)/2
st_data$sen_cl <- (st_data$sen_dem + st_data$sen_rep)/2
cor(st_data$hou_cl,st_data$sen_cl,use="pairwise.complete.obs")
pdf("party_midpoint.pdf")
plot(hou_cl~raw_mean,data=st_data,axes=F,
     xlab="State Ideal Point",ylab="State Assembly Party Midpoint")
axis(1); axis(2)
abline(lm(hou_cl~raw_mean,data=st_data))
plot(SA~hou_cl,data=st_data,axes=F,ylim=c(-.4,.3),xlim=c(-.8,.6),
     xlab="State Assembly Party Midpoint",ylab="State Assembly Indifference Point")
axis(1); axis(2)
abline(lm(SA~hou_cl,data=st_data))
plot(SS~sen_cl,data=st_data,axes=F,ylim=c(-.4,.3),xlim=c(-.8,.6),
     xlab="State Senate Party Midpoint",ylab="State Senate Indifference Point")
axis(1); axis(2)
abline(lm(SS~sen_cl,data=st_data))
dev.off()
